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MONTE CARLO EULER APPROXIMATIONS 
OF HJM TERM STRUCTURE FINANCIAL MODELS 

T. BJ6RK+, a. SZEPESSYt, R. TEMPONE§, AND G. E. ZOURARIS* 

Abstract. We present Monte Carlo-Eulor methods for a weak approximation problem related to the 
Heath- Jarrow-Morton (HJM) term structure model, based on Ito stochastic differential equations in 
infinite dimensional spaces, and prove strong and weak error convergence estimates. The weak error 
estimates are based on stochastic flows and discrete dual backward problems, and they can be used 
to identify different error contributions arising from time and maturity discretization as well as the 
classical statistical error due to finite sampling. Explicit formulas for efficient computation of sharp error 
approximation are included. Due to the structure of the HJM models considered here, the computational 
effort devoted to the error estimates is low compared to the work to compute Monte Carlo solutions 
to the HJM model. Numerical examples with known exact solution are included in order to show the 
behavior of the estimates. 



1. The HJM Model 

1.1. Generals. When valuing derivatives in the bond market it is important to use models that are 
consistent with the initial term structure observed in the market. The Heath- Jarrow-Morton (HJM) model 
for the forward rate has this property and in addition offers the freedom to choose the volatility structure, 
for example to be able to fit other derivative prices quoted in the market (see [H [71 [151 [19] ) . This HJM 
model approach is particularly suitable for Monte Carlo computations, since in general the alternative of 
tree methods leads, for the multifactor case, to non recombining trees with higher computational cost. 

In this work we focus on the numerical approximation of the price of financial instruments in the 
bond market, using the HJM model of forward rates. We propose Monte Carlo Euler methods fow which 
we develop a rigorous strong error analysis and provide rigorous weak error expansions, with leading 
error term in computable a posteriori form, offering computational reliability in the use of more com- 
plicated HJM multifactor models, where no explicit formula can be found, or such a formula is just too 
complicated to use, for the pricing of contingent claims. These weak error expansions can be used in 
adaptive algorithms to handle simultaneously different sources of error, e.g. time discretization, maturity 
discretization, and finite sampling, see [21| . To develop error estimates we use a Kolmogorov backward 
equation in an extended domain and carry out further the analysis in |21| . from general weak approxi- 
mation of Ito stochastic differential equations in R", to weak approximation of the HJM Ito stochastic 
differential equations in infinite dimensional spaces. Therefore, the main new ingredient here is to provide 
error estimates useful for adaptive refinement not only in time t but also in maturity time r. In addition, 
using the structure of the HJM model studied here, the application of a simple transformation removes 
the error caused by the representation of the initial term structure in a finite maturity partition. Finally, 
the formulas to compute sharp error approximations are simplified by exploiting the structure of the 
HJM model, reducing the work to compute such error estimates. The use of the error estimates proposed 
here is compatible with the application of variance reduction techniques, allowing for faster Monte Carlo 
computations, see [4]. 
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1.2. Description of the model. The bond market is assumed to be efficient and without friction, i.e. 
there is no arbitrage opportunity, and there exists a martingale probability measure, under which bond 
contracts can be priced as expected values of properly discounted cash flows, see [2 |31 [H]- On what 
follows, all the equations are assumed to be under such a probability measure. 

The HJM model is based on the so called forward rate, f(t,T), which relates to the price of the most 
simple type of bond, the zero coupon bond, with contracting time t and maturity time r, by 

p{t, t) = exp ^- f{t, 1]) dr] 

In particular, the non arbitrage assumption in the HJM formulation, see fT3l[T3], yields an Ito stochastic 
differential equation, for r G [0,rmax], 

d/(t,T) = Vcr^(t,T) / a^{t,s)ds] dt + y2cr^{t,T) dW'{t), tG[0,r] 
(1-1) j=i V-^* / j=i 

/(0,r)=/„(r). 

Here {W^)'-^^ are independent Wiener processes, and {a^t^T))'-^^ are stochastic processes, adapted to 
the filter structure generated by the Wiener processes. Furthermore, the initial datum for the term 
structure, : [0,Ti„ax] — >• M, is a deterministic function in C^([0, r„iax])- In this setting, the short rate, 
r{t), is defined as r{t) = f{t,t). 

On what follows the volatility function a = (ct^ , . . . , ct ') is assumed to be of the form 

where ^ : M — ^ M and A : [0,tinax] x [0,Tmax] ^ arc given bounded functions on C""''(M) and 
C™°([0,tmax] X [0,Tniax]), respectively, for toq a sufficiently large integer. Then, setting 

P = { (i, r) e [0, tniax] X [0, T,„ax] : t < t} 

problem (|l.ip reads as follows: find / = f{t, t) : V ^ M. such that 

df{t,T)^e{f{t,t))X{t,T)dt + af{t,t))X{t,Tydw{t), te[o,T], 

/(0,r)=/„(r) 
for T e [0,T,„ax], where 

(1.3) X{t,T) = Xit,T)- f X{t,z)dz, VtG[0,T], VrG[0,w]. 



Here the notation a ■ b denotes the standard inner product in R', i.e. a-b = X^^'^i a,j bj. In many models 
used in practice, the function A has the form X{t,T) — Xg{T — t), and then A(t, r) — Xg{T — t) with 

A„(T-t)^A„(T-t). / X,iz)dz. 

Jo 

Observe that to solve for / it is enough to have A^ : — M. However, in this work the usual domain of 
definition P of A and /, extends to the set [0, tmax] x [0, Tmax], leaving f\j, unchanged. The extension of T> 
helps to develop a posteriori approximations for the time and maturity discretization errors, depending 
on a linear backward problem (cf. Theorem 14. II) . 

A typical contract to price is a call option, with exercise time imax find strike price K , on a zero coupon 
bond. Its price can be written in terms of the forward rate as 

E [ e- ^0--- max {e" _ q| 

Another basic contract is a continuous cap, with price 

E / (/(t,t) ~rc)+ rfi 
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Yit)= / f{s,s)ds, Z{t)= / F{Y{s))U{J{s,s))ds, 
Jo Jo 

A{w) = / w{t) dr, W ^ L^{Ta,Tn^a.yi)- 



where Tc is a given value associated with the contract. With this motivation, and bearing in mind other 
possible contracts, wc consider the approximation of the quantity 

(1.4) E[^(/)] 
where the functional J-{f) is given by 

T{f) =f(^1 /(^: ^) ds^ G (^j^ vl,(/(t„,,,, r)) dT^ + ^ F (^^ /(s', ,s') ds'^ U{f{s, s)) ds 

with Ta being a given positive number such that < tmax < Tq < Tmax- Obviously, J-{f) is written 
equivalently as 

(1.5) Hf) = FiYit^,^)) G(A(*(/(W,-)))) + ^(tmax), 
where 



(1.6) 



The functions F : R ^ M, G : K ^ R, * : R -j- E, f7 : M ^ K, and their derivatives up to a sufficiently 
large order m^, are assumed to have a polynomial growth. We say that a function S* : R — R has a 
polynomial growth if there exist positive constants k' and C such that: |S'(a;)| < C'{1 + ) for all 
X e R. 

Let us consider the system of differential equations (jl.2p - (jl.3l) describing the dynamics for the forward 
rate / along with that for Y{t) and Z{t), i.e., 

df{t, r) = eUit, t)) A(i, r) dt + ^(/(t, t)) X{t, r)-dM^(<), 

(1.7) dY{t)=f{t,t)dt, 

dZ{t)=F{Y{t))U{f{t,t))dt, 

for t G [0,f:max] and r G [Ormax], with the initial conditions 

(1.8) /(0,r) = /„(r), r(0)=0, Z(0) = 
for r G [Ormax]- 

A approximation error for a typical discretization of the problem above will consists of a <— discretization 
error and a r— discretization error coming from the discretization of the initial condition Z^. Due to 
the special structure of (|1.7p - (|1.8p . the initial error can be avoided and practically included in the 

discretization error by introducing the anzatz 

g(t,T)=/(t,T)-/o(T), 

which implies f{t,t) = g{t,t) + fo{t). Thus, (|1.7|) - (|1.8p is formulated as follows: find g = g{t,T) : 

[0,tmax] X [0,Tmax] ^ ^ SUCh that 

dgit, t) = e{9{t, t) + /„ (t)) A(t, r) dt + (ig{t, t) + f,{t)) X{t, T)-dVK(t), Vt G [0, i„,ax], 

(1.9) dYit)={git,t)+f,it))dt, 
dZ{t)=F{Y{t))U{g{t,t)+f,{t))dt, 

for t G [0,imax], with homogeneous initial conditions 

(1.10) 3(0,r)=0, r(0)=0, Z(0) = 

for all T G [0,T,„ax]- Thus, the quantity we want to approximate takes the form 

(1.11) ¥.[g{g)] 
where 

(1-12) Q{9)--=H9 + fo)- 
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In the numerical methods, we describe later, the approximations to Y and Z will be always considered 
to be respectively the last two components of the approximate solution vector. 



1.3. Overview. Let us give an overiview of the is organized as follows. In Section [2] first we present two 
Monte Carlo Euler methods for the HJM model (|1.9|) - (|1.10|) . namely, a stochastic finite difference method, 
the Euler Finite Difference method (EFD), and a more accurate stochastic finite element method, the 
Euler Finite Element method (EFE); then, we combine a numerical quadrature rule and the outcome of 
the (EFD) or the (EFE) methods to construct a numerical approximation of the functional E [^(5)]. In 
Section[3]we provide a stong convergence analysis for the (EFD) and the (EFE) methods. Section 2] states 
and proves weak error estimates for the (EFD) method, giving explicit formulas for efficient computation 
of the discrete duals. Finally, Section [5] presents results from numerical experiments. 

2. Monte Carlo Euler Methods 

In this section first we introduce two time and maturity time discretizations of (jl.9p - (|1.10[) : the Euler- 
Finite Difference (EFD) method and the Euler-Finite Element (EFE) method. Then, we use the (EFD) 
or the (EFE) approximations along with a quadrature rule to construct approximations of the quantity 
of interest E [^(g)] defined in (|LTT|) . 

2.1. Time and maturity time discretization. Given extreme points < tmax < Ta < Tmax introduced 
in Section [I] let N and L denote the number of subintcrvals on [0, imax] and [0, T„iax], respectively. Then, 
consider partitions 

= to <•••<*« = imax and = Tq < • ■ • < = T,nax 

of the t-interval [0, tmax] and of the r-interval [0, Tmax], respectively. For technical reasons, these partitions 
are assumed to satisfy the following condition: every r-node in the interval [0,tniax] is also a t-node, i.e. 

(2.1) there exists an one-to-one index map p, such that, = t^^^j for T{ < imax- 
In addition, assume that 

(2.2) there exists an index l^, such that tmax = Te^ 
and 

(2.3) there exists an index £a such that Ta = ti^. 
Also, define the auxiliary index function, by 

(2.4) in = max {£ £ Z : <£ < L such that t/ < i„} 
introduce the notation 

At„=t„+i-f„, AWn = W{U+i)-W{t„) for n = 0,...,7V-l, 

Are = Ti+i ~Ti for £ ^ 0, . . . ,L - 1, 

and set At = maxo<n<N-i At„ and At = maxo<^<i,_i Arg. Finally, introduce the space of piecewise 
constant and right continuous functions on a r-partition, (Tg)^^^, of the interval [0,rniax], by 

S'a, = {x e L°°{0, Tmax) : thcrc are constants (q)^Zo such that xl[rf,rf+i) = ^ = 0, . . . , i - l} . 

Define the standard L^-projection 11 : L^(0,rinax) S^^ by 

/ nvxdr^ vxdr, Vx^S'^^, Vw G L^(0, Tmax), 

Jo Jo 

which satisfies 

/•re+i 

n^^|[r.T.+i) = sW / yir) dT, £^0,...,L-l, V«GL2(0,T„ax). 



For X 6 and £ ~ 0, . . . ,L — 1, denote by xe the constant value of x in [Te, Te+i). When considering a 
function, w = w{t,T), depending on two variables, the projection is always with respect to t, i.e. for 
f = 0,. . .,L - 1 and T 6 [Ti,Te+i), we have Uw{t;T) = U{w{t,-)) \ irt,Te+i) = aW /t7^' w{t,s) ds. 
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2.2. The Euler-Finite DifTerence (EFD) method. For each time level the (EFD) method approx- 
imates g{tn,.) by a piecewise constant function, g{tm .) £ S^^. In particular, it finds the approximate 
values ^„ f « g{tn, n) ioi £ ^ 0, . . . , L - 1, w Y{tn), "^^^^^^ w Z(t„) by setting first 

(2.5) %,i = ^, ^ = 0,...,L + 1, 

and, then recursively, for n = 0,...,A^ — 1, define 

fn+l/ =fn/ + ^^«C^(f«,£„ + foi'tn))Htn,Ti) 

^2g^ +?(!„/„ +/o(^«)) ^ = 0,...,L-1, 

fn+1,1. =^n,l, + At„ (f„^^^ + /o(i„)) , 

+ At„ ^^(g„^J C/(5„/„ +/o(^")) 
where the index has been defined in (12.41). 



2.3. The Euler-Finite Element (EFE) method. The (EFE) method also approximates the r-function 
g{tn, •), by a piecewise constant function g{tn, •) G S^^, but is based in a variational formulation of (|1.9p - 
(|1.10p with S^^ being the space of trial and test functions. In particular, the (EFE) is defined by the 
initial datum 

(2.7) Vo,e = 0, e = 0,...,L + l, 

and, for n = 0, . . . , iV — 1, the recursion 

9n+l,£ ~9 

+ e(f„A+/o(in)) nA(t„;r,).AW^„, € = 0,...,i-l, 

. =9n,L + [gn,l„ + /o(^n)) . 

5n+l,t + l 9n,L + l 

+ Ai„F(5„,jC/(g„^,^+/„(t„)) 
where the index £n has been defined in (|2.4p . 

2.4. Approximation of the quantity of interest E[C/(g)]. The numerical approximation of G{g) 
defined in (jl.lip involves both an approximation of the processes g, Y , Z, by computable quantities, and 
an approximation of the r- integral in (|1.6p . 

To construct an approximation of A(^((7(i,„ax, •) + /o('))) apply a composite quadrature formula, 
over the partition of [0,Tmax], based on a quadrature rule Q : C[0, 1] — R with Nq nodes Sq = {sq^i)"^-^ 
and weights Wq = (wQ.i)^j^, i.e., for v g C([0, 1]; M) the quantity Q{v) ~ X^Sl "^Q^i '^i^Q.i) approximates 
the integral v{x) dx. Also, we assume that the quadrature rule Q is of order Pq, i.e., it is exact for 
polynomials of order less or equal to pg — 1. For example, the Simpson rule has Nq = 3, Sq = (0, 1) 
and Wq = (i,|,i), with pq = 4. Another example is the Gaussian quadrature with Nq — 2, Sq = 

~ 2^3' 5 ~'~ 2^3 '^'5 ~ 5) ^^'^ P'^ ~ note that it is well known from the mathematical 

analysis of numerical quadrature that in general we have pq < 2 Nq , and the maximum value pq = 2Nq 
is achieved only by the Gaussian quadrature. 

Thus, for a fixed realization of g obtained by the (EFD) or the (EFE) method, first we approximate 
A^,{g) A(*(5(imax, •) + /o(-))) by = A(^'(£(tmax, •) + /o(0)) and then we apply the composite 

quadrature formula to construct an approximation Aq, q{'g) of A<^{g) as follows 

L-l 

A.,,Q(f)= ^ AT,Q(*(f(i 

max 1 n + ■ At£) + /„ (r£ + • An) )) 

(2.9) 



L-l 



Note that ^(tmaxj ■) is piecewise constant over the partition of [0,r,jiax] and numerical quadrature error 
in (|2.9p is caused only from the presence of the initial datum f^. In particular, if the initial datum for 
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the term structure, , is a piccewise constant function on the maturity time partition, then there is no 
quadrature error. Finally, an approximation G{g) of G{g) is computed by 

(2.10) G{V) = F{V.,,)G{KAV))+V.,.+^- 

The Monte Carlo method, [16], approximates the expectation of a given random variable X by a sample 
average of M independent realizations of X, i.e. E[X] « A{M; X) = jj X^jLi ^i'^j)- In particular, here 
we approximate E[Cy((7)] by a sample average of 5 (5), 



(2.11) A{M-g{g)) ^ ±^Y.[F{g,^,{u,)) G{K.,M^,))] 

Therefore, the exact computational weak error 

(2.12) S, = ng{g)]-A{M-M)) 
naturally separates into three error contributions as follows: 



(2.13) 

with 

(2.14) 



— Eo + Eq 



E, 



E,= nQ{9)]-^[Q{9)]. E^=V.[g{g)] -E[g(g)], 

Es = ^[gm]-A{M-g{^)) 



where E^ is the error contribution from t- and r- discretization, Eq is the quadrature error in (12. 9p . and 
Eq is the statistical error. 



3. Strong Convergence 

To carry out an error analysis for the numerical methods proposed in Section [51 we assume that there 
exists nonnegative constants Cj_i and C^ 2 such that 

(3.1) 1^2(^)1 < + V:r gK, 
and 

(3.2) \s^'^(^nc)-s:^{z)\ + \^{^x)~i{z)\< C^a\^-zl Vx, zeM. 

3.1. Bounds for moments. In Lemmas 13.11 and |3 . 2[ we show, respectively, boundness for the moments 
of the T— derivatives of the solution g to the problem (|1.9p - (|1.10p . and for the functional value g{g)- 

Lemma 3.1. Let D^, = [0,<inax] x [0,Tinax], g be the solution of (|1.9p - (jl.lOp and v E Nq. Also, we 
assume that the derivatives d^X and {d^Xj)j^i are well defined and continuous on Z?*, for £ = 0, . . . ,1^. 
Then, for £ = 0,. .. ,1/ and k S N, there exists a positive constant C"^, depending on k, £, {d^Xj)j^i, 
9rX, fa, Q,!, Tmax and fmax, such that 



(3.3) 

where C^,i is the constant in (j3.1 



max 



\dig{t,T)\ 



Proof. Let k G N, £ € {0, ■ ■ ■ ,v} and (t, r) S D*. Also, in order to simplify the notation, we set := fmax 
and T,nax- Our first step is to use (|1.9p to get 



(3.4) 
where 



\dU{t,r)r < (J+1)^ 



'TLit,r) + TUt,r)[ 



T(jt,T)^E 



TUt,r)^Y.^ 



diX{s,T)e{g{s,s) + f,{s))ds 



2k 



' diX,is,r)agis,s) + f„{s))dWHs) 



Using (j3.1l) and applying the Holder inequality we have 

' \di\{s,T)\ (l + |/„(,s)| + |,9(s,s)|) ds 

' \d[\{s,T)\{l + \fM\)ds 
+ (^j\diX{s,T)\ \g{s,s)\ds 



2k 



(3.5) 



where Cf = 22«- 



max (/„**|a^A(.,r)|(l + |/„(.)|) 



2k 

rfs ) and 



= 2^^^-^ max / |9;A(s, t)| — ds 

r G [0,r ^ 



2k-1 



Next, using the properties of the Ito integral and p.ip . we obtain 

TUt,r)<i2K-l)\lj2([\d%is,r)fE[ei9{s,s) + Us))] ds 



(3.6) 



< 



< 



{2k ~ 1)!! (Q,i)« ^ f / (9^A, (s, t)Y (1 + |/„(s)| + E [ |.9(,s, ,s)| ] ) ds 
j=i 

(2«:-1)!!(Q.i)«E(/(^'^^(^'^))' (2+|/„(s)|+E[|.g(s,s)n) ds 
l[\g{s,s)\''] ds 



where C3-'^^ (2k - 1)!! 2«-i (Q,i)'^ (e;=i max |af A, ^ j and C^^ ^ C^'' (2+|/„(s)|) 
Now, combine p.4p . (|3.5p and p.6p . to arrive at 



ds 



(3.7) 



E 



{d'Mt.rf^' <C^'' + C^/[ I E[|5(s,s)n ds) + C^;/ I E[|g(s,s)r] ds 



where Cf = ( J + l)2--i ( + ), = ( J + l)2«-i and C^/ = ( J + l)2--i 
Consider the case k = 1 and ^ = 0, and set r = t in p.7p . to obtain 

E[|g(i,t)n < C^" + (C;;" + Ci;°) T E [ |5(s, s)^ ds, Vie[0,i.], 

Jo 

which, after the application of the Gronwall lemma, yields 

(3.8) E[|.g(t,t)|2] < C>Oe(^i"+^i")*, Vte[0,t.]. 
Now, combine (|3.8p and p. 71) (with k = 1), to get 

(3.9) E\\dig{t,r)n < + '^^^^M^ Uc\f'^c]-,t _ ,1 ^(,^^) ^ 

L ' 'J '-'II +'-'111 L J 

for ^ = 0, . . . , I', which establishes p.3p for k = 1. 

Now, consider the case k > 2. Then, use p.Sp and p.7p . to obtain 

(3.10) E[|9;5(t,T)f"l < C;^^ + C«/ /*E[|.g(s,s)p«] ds, V(i,T)ei?., £ = 0,...,i., 
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where = ' + 
obtain 



■^11 ' ^ III 



Take £ = and set r = i in ([XTU)) . to 



E[|g(t,i)l'"] < + / E[|g(s,s)|2'^] ds, Vie[0,t,]. 

Jo 

Apply again the Gronwall lemma, to conclude that 

(3.11) E[|5(t,i)|"'] < C;^°eC'™ *, Vie[0,i.]. 

Finally, combine p.lip and (|3.10l) to have 



(e^™*-l), V(t,r)£i^,, £ = 0,, 



which yields the desired bound p.3p for k > 2. 



□ 



Lemma 3.2. Let {g,Y,Z) be the solution of the system (|1.9p - ()1.10p . ylfeo, we assume that the functions 
-F, G, [/ : M — > M have polynomial growth pp, pa, p^ and pu with constants Cp, Cq, and Cu, 
respectively. Then, for k £ N, there exists a positive constant C^, depending on k and the data of the 
problem, such that 



(3.12) 



|i^(i"(tmax))r 



-E 



|G(A(*(g(W,-) + ./o)))r 



E 



l^(imax) I 



Proof. Let K G N. To simplify the notation, we set := Tmax, := ^max and T(t) ■= g{tnia,x,T) + f^ir) = 
fitmaxiT) for T £ [0,Tmax]. Sincc F, U and G have polynomial growth, using the Holder inequality and 
(|3.3p for £ = 0, we obtain 

E[|F(r(t))|2'^] < (C^)2« 2^'^-! (l+E[|y(t)|2«PF]), VtG[0,t,], 

E[|y(i)|2"] <(2i.)''"-i /*(E[|5(s,s)|2"]+|/„(s)|2") 

Jo 

< (2i.)2"-i / * (C|;„,o + l/o(s)l"") rfs, Vt G [0,t,], Vm G N, 



(3.13) 
(3.14) 

(3.15) 
(3.16) 

(3.17) 
(3.18) 



E[\FiY{s))\'^] ds+ / E[|C/(/(s,s))r] , 



E[|Z(i,)|2«] <it,f''- 
E[|C/(/(t,0)P"] <(a)2"32'"-i (l + |/„(t)p™Pc. +E[|g(f,t)p'"Pc.]) 

< (a)2" 32™-i ( 1 + |/o(t)P'"P- + G2%^^,o) , Vi G [0,tj, Vm G N, 
E[|G(A(*(T)))|2«] < (Caf- (l+E[|A(*(T))|2-fG]), 



E[|A(vI/(T))|2'"] <{r^-raf 



E [|*(T(t))|2"1 dr, VmGN, 



and 
(3.19) 



E[|vI/(T(t))|2"] <(G,)2'"32"-i (l+E[|5(t.,r)|2™P*] +|/„(t)|2"p*) 

< (G,)2" 32—1 ( 1 + |/,(r)|2"f* + G2%^^„o ) , Vr G [0, r.], VmGN. 
Thus, we obtain p.l2p combining the inequalities p.l3p - p.l9p above. 



□ 



In Lemma 13.31 below, we show boundncss for the moments of the numerical approximations produced 
by the (EFD) and the (EFM) method. 

Lemma 3.3. Let I :~ {0, . . . , N} x {0, . . . , L — 1} and (g„ i)(n.i)ei ^^'^ numerical approximations 
produced by the (EFD) or the (EFM) method. Then, for k G N, there exists a nonnegative constant C^\, 



depending on k, {Xj)j^i, A, /„, G^^i, 
(3.20) 

where G^_i is the constant in (|3.1I) . 



and imaxj such that 



max E n 



2k " 



< G" 



Proof. Let = [0,imax] X [0,rmax], K e N, {n,£) e I with n > 1. Then, from (US]) and (E^, we 
conckide that 



(3.21) 



m=0 



m— J — 1 



where /™ /„(t,„), i."^^ = A(t™, r^) and /i^'^ = Aj(i„, ) for the (EFD) method, and t^"^^ = nA(t™; r^) 



and /i"''^ = IlXj{t„i;Te) for the (EFE) method. Thus, we obtain 



(3.22) 
where 



2k 



2.K E' 



=0 / 

V m=0 / 



2k;' 



Using p.ip wc bound T^'^ as foUows 



5] Ai„,|^."'^l(i + l/;^l + l5,„.,„|) 



v,m=0 



2k 



2k 



Km=0 



which, after applying the Holder inequality, yields 



(3.23) 



\9m,e„ 



|2k 



where Cd.i,k = (2 C^.i)^" (tmax)^''""'^ max^,^ [ |A| (1 + j/^ |) J^'^. Also, using the properties of independent 
Gaussian random variables and p.ip , we obtain 

.7 / n-l \« 

Tr-^<(2«-l)!!5] ^At™(^7-^)2E[e^(f„,,„+/™)] 



<(2«-l)!!(Q,i)«m^ax|Ap'^ ( ^ At™ ( 2 + + E [ ^ ) 



ri-l 



, m=0 



<(2k- l)!!(Q,i)'' maxlAp'^ 
which yields that 



Tl-l 



tn.ax max (2 + l/J) + ^ At™ E [ |g™ 



(3.24) 



±2 S ^D,2,K 



n-l 



(imax)'^+ ^ At™E[|g™_,^|2] 



, m=0 



where Cn,2,K = 2"^^ {2k - 1)!! (Q,!)*" max^^ |A|2« max[o,f.„^^j (2 + l/J)". Now, combining ([3^, ((3:231) 
and p.24p we obtain 

n-l / n-l \ 

(3.25) E [ |f„,,p« ] < C;-- + C-'- ^ At™ E [ |f™,,^ ^ ] + C^'" E [ |f™^,,„ , 

ni=0 \ m=0 / 

where C^ °, C" " and C™ " are constants that depend on J, k, tmax, Cd,i,k and Cd,2,k- 
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First, let us consider the case k = 1. Then, setting £ = i!„ in p.25p . we obtain 

n-l 

(3.26) E[|f„_,j2] < C(-°+Cr'" ^ At™E[|f,„_,^j2], n = l,...,7V, 

m=0 

where Cf ^ = Cf ^ + C"'-"". Setting ;3„ := E [ |f „ ^ pi for n = 0, . . . , iV, is written equiva- 

Icntly as follows 

n-l 

(3.27) /3„ < 1 + Cf'" ^ At,n Pm, n = l,...,N. 

Now, setting pi := 1 and p„ := 1 + Cf' " X]m=\ ^^'n Pm ioi n — 2, . . . , N and observing that Pq — 0, we 
use (|3.27p and apply a simple induction argument to get 

(3.28) /3„<p„, n = l,...,N. 

Since p„ = (1 + C(^'° At„_i) p„_i for n = 2, . . . , A^, wc use the inequality > 1 + x for x > 0, and a 
simple induction argument to conclude that 

(3.29) pn < cxp(C7r" i„), n = 1, . . . , Af. 
Thus, ((X^ and (j^:^ yield 

(3.30) max E [ |f„ ^ ] < exp(Cr" U,^), 

0<m<N m J 

which, along with p.25p . establishes p.20p for k = 1. 

Now, we assume that k > 2. Then, we combine p.25p and (|3.30p to obtain 

ri-l 

(3.31) E[|f„^,j2«] < ^ At„,E[|f,„_,^j2«] , n = l,...,7V, 

where C^-" = +C™'" (^max)" (Cf"' exp(Cf ^ imax))''. Then, proceeding as in obtaining ([5301) from 
p.26p . we arrive at 

(3.32) max E [ , ] < C^" cxp(Cf i^^x), 

0<m<iv 

which, along with (P?^ and (13301), yields ((X^ for k > 2. □ 



3.2. Estimates for the consistency error. In Lcmmas l3.4l and l375l bclow. wc show that some Lipschitz- 
typc properties for the solution g to the problem (|1.9p - (jl.l0p hold. 

Lemma 3.4. Let k G N and g be the solution of (|1.9p ~ (|1.10p . Then, it holds that 

(3.33) E[|.g(f,ri)-g(t,T2)P'=] < C^^i|ri-T2|2«, Vn, r2 6 [0, w], Vt £ [0, W], 

where C"i is the constant in p.3p for £ — I. 

Proof. Let t G [0,tniax] and ri, T2 G [0,Tmax] with T2 > ti. Then, applying the Holder inequality, we have 



E[\g{t,n)-g(t,T2)\^^] =E 



, 2k 

drgit, t) dr 



Tl 



<|t-2-ti|^'^-^ / E[\drg{t,T)\'-\ dr 
<|t2-tiP'' max E [|a^g(t, r)p'' ] . 

re[Ti,T2] 

Thus, wc obtain p.33p combining the inequality above and (|3.3p for £ = 1. □ 

Lemma 3.5. Let k G N and g be the solution of (|1.9p - (jl.l0p . Then, there exists a nonnegative constant 
Clip, depending and tniax7 suc/i t/iat 

(3.34) E[|5(ii,T)-g(i2,T)|2«] < CL.p|ti-t2r, Vti, t2 e [0, tn,ax], Vtg[0,w]. 
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Proof. Let = [0,iniax] X [0,rmax], T € [0, T^ax] and ti, ^2 e [0,tmax] with t2>ti. Proceeding as in the 
proof of Lemma 13.11 we obtain 

(3.35) E[|5(T,ti)-5(T,<2)|'"] < [B,,,{T;h,t2) + B^^j,{T;h,t2)] 

where 

rt2 



A(s,T)e'(5(s,s) + /„(s))rfs 
rt2 



2k 



Using the Holder inequahty and (|3.ip we obtain 



S„,,(r;ti,t2) <(2Q,ir'' maxlAI^'^E 



(3.36) 



and 



(3.37) 



t2 



{l + \Us)\)ds + \gis,s)\ds 



<(2Q,i)2- maxlAl^lti-ial'" max (1 + |/„(s)|)2- + max E [ |5(s, s)|2« ] 

\se[tl,t2] sG[tl,t2] 



5«,z.(r;ii,t2) <(2k- l)!!(Q,i)« ^niax|Aj 



2k 



t2 



(2+|/„(s)|+E[|5(5,5)n) 



< (2«-l)!!(Q,i)'= K]max|A,f« \t2 - max (2+\f,{s)\+E[\g{s,s)\'])\ 

1 — ' se|ti,t2j 



Thus, ((334)1 follows easily from (jXSSl) . ([3361) . ([XST)) and ([331) for f = 0. 



□ 



In Proposition l3.1l that follows, we prove a consistency result for the (EFD) and (EFE) methods defined 
in Section [2] 

Proposition 3.1. Let k e N, g be the solution of ([Q)) - (jl.lOp . f^ := f„{tm) for m = 0, . . . , iV, and 

K,n,i be defined by 

(3.38) gn+i,i = gn,t + A<„ z."^^ e'(?n,c + /;) + 51 m"'' ?(5n,^„ + /:) Aiy^ + /C„,,, 

/or n = 0, . . . , iV — 1 and ^ = 0,...,L — 1, where z^"'^ = A(t„, r^), ^J'^ = Xj{tfn, t^) and gm.e = g{tm, T() for 
the (EFD) method, and i/™^^ — IlX{tm;Ti) and fi^'^ = llXj{tm',Te) and gm.e = ^g{tm',Ti) for the (EFE) 
method. Also, we assume that f^ € C^([0, Tmax]; K) and dtX, dtrX, {dtXj)j^i, {dtrXjY^^i are well-defined 
and continuous on [0,tinax] x [0,T,„ax]- Then, there exists a nonnegative constant Ccn,i, independent of 
the partitions of the intervals [0,tinax] o,nd [0, Tmax]; such that 



(3.39) 



J2 ^r, 

m=0 



< CcnA [(Ai)" + (Ar)2-] 



for n = 0, . . . , N ^ I and £ ~ 0, . . . , L — 1. In addition, for the (EFD) method there exists a nonnegative 
constant Ccn.2, independent of the partitions of the intervals [0,iinax] o,nd [0, Tmax]? such that 



(3.40) 



E 



m=0 



< Ccn.2{^nf^ [(Ai)'= + (Ar)2-] 



forn = Q,...,N-l and £ ^ 0, . . . , L - 2. 



Proof. Here, wc set D^, := [0,imax] x [0,T„iax] and use the symbol C for a generic constant independent 
of the partitions of the intervals [0,tinax] and [0,rinax]- First, we observe that (|1.9p yields that 



(3.41) 



e=9n.e+ / v''is)^''{g{s,s) + f,{s)) ds 



1=1 



forn = 0, ...,A^ - 1 and £ = 0,...,L - 1, where v^{s) = \{s,ti) and /^^(s) = \j{s,n) for the (EFD) 
method and v^{s) = I{\{s;ti) and /xj(s) = nAj(s;rf) for the (EFE) method. Then, subtracting p.4ip 
from p.38p we obtain X]m=o ^rn,i = Ylt=i ^?,c for ^ = 0, . . . , L — 1 and n = 0, . . . , N — 1, where 



m=0 



m=0 
j=l m=0 
j=l m=0"^*'" 

Next, using p.ip . the Holder inequality and (|3.3p . we obtain 



{l + \f„{-'^)\ + \g{s,s)\)ds 



(3.42) 



<C(At)2 



(l + |/„(s)|ds + 



2k 



\g{s, s) \ ds 



2k 



<C(Ai)2« 

<C(At)2« [l + (Uax)'"C^:o 



l + (tmax)^'' max E[(5(s,s)) 

Se[0,tmax] 



2k " 



and 



E 



iE^;^r] <CE 



E 1 {l^]{s)~y^';''')(,{g{s,s) + m)dW^{s) 



(3.43) 



J — 1 \_rn—Q 

2k 



^^E 

< C (At) 
<C(At)2'' 
<C(At)2« 



T?J=0 *" 



E W - M?')' IE [C'(.9(s,s) + /„(«))] 



(l + |/„(s)|+E[|,g(s,s)|]) ds 
(2 + |/„(s)|+E[(g(s,s))2]) 
(2 + |/„(s)|+C7ro ) ds 



Now, we apply p.2p and the Holder inequality, to get 

" /•t„^+l 



E 



<CE 



<CE 



(3.44) 



E / {\9i^,^)~9m,iJ + \fois)-fr\) ds] 



r,Jt 



^ m— 



2k 



^ 771 — 



<c 



<c 



m=0 •^*''» 



(Aif« max ^ / E[|g(s,s)-g™.,„r] ds 



and 



E 



" ft 



m—O 



E / A^?' [^(9is, s) + f,{s)) - a97n,i„. + /r)] dW^s) 



2k 



(3.45) 



<C 



<C 



<c 



E/ (A^r'fE (e(3(s,s)+/cW)-^(W„+/r))^ 



m=Q 



E / (E[|5(s,s)-ff™,,„j2]+|/„(s)-/„"p) ds 
(At)2 max |/„f + E / E[\g{s,s)-g7n,ij'] 



Using ([OS]) . ([3:341 and we have 

E [ \g{s, s) - ] < C (e [ |g(s, s) - ^(s, t,,JP« ] + E [ \g{s, nj - .9(t^, ] 

+ E[|.g(i,„,T,„) -5„,,,,J2-] ) 

(3.46) < C ( |S - -^C. 1^" + |S - ^mT + E [ |.9(<m, Tf„J - 5m,£„. I^'' ] ) 

< C (|S - t™|2K + 1^^^^ _ |2k ^ (^^)K ^ ^ [ |^^(^^^ ^^^j _ ] ) 

< - I'" + (At)" + E [ |5(t„„ nj ~ |2« ] ) 

< C ((Ar)2- + (At)« + E [ |,9(<„, ) - P'' ] ) 

for s G [tm, im+i] and m = 0, . . . , — 1. For the (EFE) method, after using p.33p . we have 



E[l5(*m,T-f„)-5„,f„p-] <(Arf„)-i / " E[(g(i„,r,„J-g(i„,r))2«] dr 



(3.47) 



<C(Ar,„J-i / |r,„_r|^-dT 

<C(Ar)2«, 77i = 0,...,iV-l, 
while for the (EFD) method the term we estimate above vanishes. Finally, p.46p and p. 471) yield 
(3.48) E [ |.g(s, s) - ] < C [(At)'' + (At)^^] 

for s e [t„„ t„,+i] and m = 0, . . . , iV - 1. Observing that E \{YJ'i=o ^m,if''\ < ^^"'^ Eti ^ [(^^"c 
for ^ = 0, . . . , L — 1 and n = 0, . . . , — 1, and that estimate p.48p holds for k = 1, the estimate (|3.39p 
for the consistency error follows easily in view of (P^i^ . (|3ZS1), dSilll), (P^^ and 
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Since, E [(ELo " ^-m+i) J < 4^'^-^ Eti E [{El;^^' - iM;'?'' 

= 0,...,^-l, we obtain ((5^ for the (EFD) nethod, observing that 



for ^ = 0, . . . , L - 2 and 



|ym,f+l _ym/| <CArf 



(3.49) 



where y^' 



<CAtATe, Vs e 



or /i"'^ and y = A or Xj, respectively, and proceeding as above. 



□ 



3.3. Error estimation. In this section we derive an error estimate for the strong approximation error 
Q{g) — Qig) by spHtting it as sum of the strong discretization error G{g) — Q{g) which we estimate in 
Theorem [32] and of the strong mumcrical quadrature error Q(£)~Q(£) which we estimate in Theorem l3.3l 



Theorem 3.2. Let g, Y and Z be the solution of P3|) - P^ . := {0, . . . , TV} x {0, . . . , i - 2} 
J := {0, . . . , N} X {0, . . . ,L — 1}, I := {0, . . . , N} x {0, . . . , i + 1} and {£n,i)(n,t)ei be the numerical 
approximations produced by the (EFD) or the (EFM) method. Also, we assume that the functions F, 
F' , G, G' , U, [/' : IR ^ R have polynomial growth, and we define Aq,{w) :~ A{'i'{w + fa)) for w £ S^^ 
or w € C([0, Tijiax]; 1^)- Then, there exist nonnegative constants (Cf^ )^^^, independent of the partitions 
of the intervals [0,ti„ax] a'^'^ [Oi'^'maxlj such that 



max I E 



0<n<N 



(3.50) 
(3.51) 
(3.52) 
(3.53) 
(3.54) 

and, for the (EFD) method, 



\gitn,Ti) ~ g^f?'^ 



~ < [(At)5 +Ar], 

max ^ (E [ |r(t„) - f r ] ) ^ < Cf [ (At) U Ar ] , 

E[l^(W)-f„,.+,r])* < C3-'[(At)UAr], 

^' < Cr[(At)^+AT], 



E 



\A^,{g)-A^{g) 



{E[\gig)~gm'-])^'- < Cr[(At)UAr] 



(3.55) 



max I E 

{n,e)eM 



< crim^+Ar]. 



Proof. Here, we set := [0, tmax] x [0, Tmax] and wiU use the symbol C for a generic constant independent 
of the partitions of the intervals [0,tniax] and [0,Tniax]- Let Em,i = dm,i — Vm ^ for m = 0, . . . , and 
£ = 0,...,!/ — 1, where g„i^i = g{tm,Te) for the (EFD) method and g^^e = ^g{tm\Tg) for the (EFE) 
method. First, subtract (|2.6p or (|2.8p from p.38p . and then sum with respect to n, to obtain 



En,i = An,e + Bn,e + ^ }Cm,e, £ = 0, . . . ,L - 1, n = l 



,...,iV, 



m=0 



where 



n-1 

m— 

71—1 J 

m— j — 1 
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/r - /o(*™)' = Htm,Te) and m™^' = A,(t™,r,) for the (EFD) method and 



(3.56) 



UXj{tm',Ti) for the (EFE) method. Thus, we have 



E [ |£;„,,|2- ] = S^-^-i E [ |A„,,|2« ] + E [ ] + E 



(3.57) 



and 



(3.58) 



,m,e 



2k 



Ii\{t„T,;Ti) and 



for ^ = 0, . . . , L — 1 and n = 1, . . . , iV. First, using p.2p and the Holder inequaUty, we obtain 



2k 



ri-1 



<C ^ At,„E[(i;^,,„)2«] 



ri-l 



J — 1 \ m— 
/ n-1 



, m=0 

' Tl-1 



<C ^ At^E [(£;„,,„ )2] 



for £ = 0, . . . , L - 1 and n = 1, . . . , iV. Combining, (|336)) . ((3:57)) and (|338)) and (|339)) . we have 



(3.59) E[|£;„,f|2«] <C 



((At)« + (At)2-) + ^ At„, E [ ] + At™ E [ p ] 



?n— 



^Tn=0 



for i = 0, . . . , i — 1 and ri = 1, . . . , A^. Considering the case k = 1 and proceeding as in the proof of 
Lemma [3731 from p.59p we arrive at the estimate 



(3.60) 



max E\\En,e„f] < C{At+{ATf). 

0<n<iv 



Letting k > 2, under the view of p.60p . the inequahty p.59p yields 



(3.61) 



n-l 



{{Atr + (At)2«) + ^ At,„ E [ ] 



for £ = 0, . . . , L — 1 and n = 1,. . . ,N. Now, proceeding again as in the proof of Lemma 13. 3( from p.6ip 
we conclude that 



(3.62) 



max E\\E„e l'^^] < C ((At)'= + (At)^'') 

0<n<JV L I . n I J \v / \ / / 



Thus, combining p.6ip and p.62p we arrive at 
(3.63) 



max E [ |£;„.^P'=] < C {{AtY + (At)^^) 

n,f)eIjv,L 



The estimate (|330| for the (EFD) method follows directly from (|3:63| . For the (EFE) method, (|330| 
follows combining p.63p with the following estimate (cf. p.47p ) 

max E [ |5(t„,rf) -ng(i„;r£)l2''] < max max E [ |5(t„, r) - g(t„, r^)!^'' 1 



-A„ I . Br 



<C(Ar) 



,2k 



for ^ = 0, . . . , i - 2 and n = 1, . . . , A^, 



Since ~ Ati ~^ Are L^m=0 Art „.r-r^ 

to obtain the estimate p.55p for the (EFD) method we proceed as above using p.40p and p.49p . 
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In order to get the second error estimate, we use (|1.6p and (j2.6|) or (j2.8p . to conclude that 



(3.64) 
where 



E 



-1 



-1 /.t 



n-l 



m=0 



First, we observe that 
(3.65) 

Next, we use the Holder inequality, p.34p . p.33p and p.50p to obtain 



|Cri^''< CiAtr max^|/J«, n = l,...,7V. 



^1=0"^* 



< 



^E 



-1 pt. 



[\gis, s) - g{t^,Te Jl^"] ds 

[|,g(s,s) -5(t,„,s)P« + \git,n,s) - giUn^nJl^""] ds 



(3.66) 



<C 



<C 



(At)^ + E At™ |r^„+i 



Tin 



m=0 



and 



(3.67) 



<C [iAt)^ + {ATf^], n^l,...,N, 



E[|C3"l'1 <C AUr^E[\g{t„„T,J ^g^jj'^] 



<C [(Ai)" + (AT)2-] , n = l,...,7V. 

Thus, follows easily from ^M^, <^Ml, and dSU]). 

In order to prove our third error estimate, we use (jl.6p . (|2.6p or (|2.8p . and the mean value theorem 
for scalar fields, to conclude that 



(3.68) 
where 



^(tmax)-f™r1 < CElE[|r,p-], 



pt^ + l 

Ti E / ^'(^™ f/(S™(s)) iY{s) - g^J ds, 

m=0 

/•*,„+! 



^2 E / F(A™(5))C/'(B™(s))(g(s,s) -g„_,^) ds, 
Ta E / C/'(S™(.s)) (/o(.) - C) ds, 



and 



Amis) :=(5™(s) (r(s) - .g„^^) + 

Sm(s) ■.= 6m{s) ig{s,s) + f„{s)) + (1 - ~5^{s)) (g^^,^ + /™), 
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with 6m{s); Sm{s) £ [0,1]. Let e N. Since F, F', U and U' have polynomial growth, we use p.Sip . 
p.l2p . (|3.20p and p.3p to conclude that there exist a nonncgative constant C™ such that 

(3.69) max sup [E [ |F'(A„,(s)) ;7(B™(s))|2™ ] + E [ |^^(A,„(s)) ] ] < Cf . 

0<™<"-ise(t„,t„+i) 

Also, we use the Holder inequality and p.3p to arrive at 

E[|y(t,)-y(tjp*] <(t,-ij2n-.-i rE[|.g(.,5) + /„(.)p*] ds 

(0.7U) Jta 

<C(4-ia)'™ 

for all ta, tf, S [0,iinax] with ta <t}j. Now, we are ready to estimare the quantities at the right hand side 
of p.68p . First, we use the Holder inequality and p.69p to arrive at 

E[|rip-] <C / E[|^^'(A„,(s))C/(B„,(s))p«|r(s)-f,„^j2K] 
m=o 

(3.71) <C Xi r'"^'(E[|^^'(A„(s))[/(B™(s))|4-])^ (E[|y(s)-f„^,r])^ ds, 

m=0 

<cE / (iE[|r(s)-r(ur + |r(t„)-f„,j'^''])^ ds, 

,71=0 

E[|r2p«] <C r"^\E[|F(A„(s))C/'(i?,„(s))r«])^ (E[|5(5,s)-f„,,^r])^ ds 



(3.72) 



m=0 ■ 

(-tm+l 



<^E/ (E[l5(5,5)-5(t™,T,„)|4« + |5(t™,T,„)-5™,£„l'1)" c^s, 

m=0 



and 

(3.73) E[|r3p«] < C r^'\.fo{s)-fr\'^ds. 

m=0 

Next, we combining (|3J2l) . ((3J3)) . (1X48)) and ((3?50)) we obtain 

(3.74) E[|r2|^''] + E[|r3|2''] < C [(At)'' + (At)2«] . 
Finally, we combine p.7ip . p.70p and p.5ip to obtain 

(3.75) E[|ri|2''] < C [{At)^ + {ATf''] . 

Thus, the error estimate (|3.52p is a simple consequence of (|3.68l) . p.74p and p. 751) . 

To derive our fourth error estimate, first we set E :~ A^{g) — Aq,(g), and then wc use the Holder 
inequality to obtain 

(3.76) E[|^|2-] < C/bI/b;, 
where 



sup l^-' (/„(r) + eg(W,r) + (1 - e)ff„,^)|'*" 



,££[0,1] 

E.-.^Y. / E[|5(W,r)-f„,,|4''] dT. 

Since has polynomial growth, the use of p.20p and p.3p yields that 

(3.77) Ea < C. 
Also, using p.33p and p.50p we obtain 

(3.78) E^ < C [ {Atf^ + {Arf'^ ] . 

Thus, the estimate p.53p follows after combining p.76p . p.77p and p.78p . 
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dr. 



To obtain our fifth error estimate, first we set Eg Q{g) — Q{g) and use the error bound (|3.52|) to 
obtain 



(3.79) 
where 



g,, :=E[|G(A,(g))| 



g., :=E[|i^(y(Wx))-n5«,Jl'1 
g^, :=E[|G(A„(.g))-G(A,(f))ri. 



Since F and G have polynomial growth, we combine p.l2p and p.5ip to get 

(3.80) g^^ + g^, < c. 

Since F' has polynomial growth, we use the mean value theorem, the Cauchy-Schwarz inequality, (|3.12p 
and the error bound p.5ip to have 



(3.81) 



^E 


max 




.ee[o,i] 



i^'(er(W) + (l-e)ff. 



I 8k 



(E[|y(w)-5«,J'1)^ 



< 



C [(Ai)2« + (At)^'^] 



Similarly, since G' has polynomial growth, we use the mean value theorem, the Cauchy-Schwarz inequality, 
p.l2p . and the error bound (j3.53p to have 



(3.82) 



^E 


max 




.ee[o,i] 



\G' (eA,(.9) + (1 - 6) A,(f)) I'^^J j (E [ |A,(5) - A,(f)|«'^])^ 
Thus, the error bound ([334]) is a simple consequence of (|3?79l) . (|3^ . ([3?8T|) and (|3^ . 



□ 



Theorem 3.3. Let T :~ {0,...,A^} x {0, . . . , L + 1}, {gni)(n.i)ei ^'^ numerical approximations 
produced by the (EFD) or the (EFM) method, A^{'g) be defined as in the Theorem \3.S\ and Q{'g) 
be the quantity defined by (|2.9p . We assume that the quadrature rule Q used in (j2.9p is of order p^ , 
* G C^Q(M;R) and /„ e CPq([0 J ''"max] i IK.) . Also, wc assumc that ^ and all its derivatives up to order 
Pq, along with the functions F and G' , have polynomial growth. Then, for k G N, there exist constants 
C^'^ and G^'^, independent of the partitions of the intervals [0, tmax] o,nd [0, Tmax]; such that 



(3.83) 

and 

(3.84) 



(E[|A,(f)-A,,«(f)f^])'^ < Gr(Ar)^Q 
m)-mT]T < (Ar)^«- 



Proof. For ^ = . . . , i — 1, we set vg{s) := ^'(^jy ( + f„{Te + s At^)) for s G [0, 1]. Since the quadrature 
rule Q has order p^ , applying a standard argument from the error analysis for quadrature rules based on 
the Taylor formula (see, e.g., [5]), we obtain 



(3.85) 



E 



A45)-A,,^(g)r^ < G(At) 



2k, 



E 



max max 

ea<i<L-i [04] 



2k 



Observing that ds'' Viis) = ^'^o "^^'HV N.e + foin + s An)) /r« "(n + sAn)), assuming that 
has polynomial growth pj for j = 1, . . . and using p.20p . we obtain 



(3.86) 



E 



max sup 

,<£<i-i [0,1] 



<G. 



Now, combine ((5^ and ((5^ to arrive at ((5^ . 
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Since F and G' have polynomial growth, using the Cauchy-Schwarz inequality, the mean value theorem, 
(ESD), (EHI), and (jXS^ . we obtain 



E 



<C E 



\F{9...)\ 



E 



<C IE 



sup |G" (eA^(g) + (1 - e) A^,q(5)) 
£e[o,i] 



|A*(5) - A*,q(<?)|*^'' 



<C(Ar)2'=PQ 
which yields the estimate p.84p . 

4. Computable Weak Error Approximation 



□ 



In this section we present a computable approximation for the weak t— and r— discretization error 
Eo defined in (|2.14p for the (EFD) method. In Theorem 14.11 below we give an estimate of Eo which, as 
the step size of both the time and maturity time partitions go to zero and the number of realizations 
goes to infinity, is asymptotically correct. On the other hand, the statistical error Eg can be analyzed 
by the Central Limit Theorem or Berry-Esseen Theorem, a standard procedure in Monte Carlo methods 
(cf. Scction[S]). While, in Theorcm l3.3l we have estimated the quadrature error Eq, concluding that when 
the order Pq of the quadrature rule Q we use in (|2.9p is sufficiently large, the quadrature error, Eq, is a 
higher order term in the expansion of the computational error. 

To have an easier access to the results and the techniques of [H], we reformulate problem (|1.9p - (|1.10p . 
letting the process g — g{t, r) be the solution of the problem 

dg{t,T)^a{t,T,g{t,t))dt + b{t,T,g{t,t))-dW{t), Vt £ [0, tn^ax], 

5(0,t)=0, 

for T e [0,Ti„ax], where a : [0,t,„ax] x [0,T„iax] x M ^ R, 6 ; [0,ti„ax] x [0,r„iax] x M ^ M' given by 

a{t,T,x) =e{x + f„{t))\{t,T), 

b{t,T,x)=aX + f,{t))X{t,T). 

We approximate the unknown process g{t,T) by a time and maturity discretization g{t,T), with t e 
{tn)n=o ^'^d r e {Te)iZQ, based on the (EFD) method, which, for n = 0, . . . , — 1, reads 

g{tn+i,n) =g{tn, Ti) + a{tn,Te,gitn,Te„)) Atn + b{tn,Te,g{tn,nJ)-AWn, e = 0, . . . ,L ~ 1, 

f(0,T£)-0, ^ = 0,...,L. 

For the analysis of the (EFD) method, it is useful to extend its definition for all times t and all maturities 
r as follows: for n = 0, . . . , A^ — 1 and £ = 0, . . . , L — 1, set 

f(t,r) ^^{tn,Ti) + a{tn,Ti,^{tn,Ti,J){t - t„) + 6(t„, T£, f (t„, T£„ )) • {W (t) ~ W{tn)) 

(4.3) =^{tn,Ti) + W{s,T,'g{tn,TiJ)ds + S(s, T, f (i„ , T£„ ) ) •dW^(s) , V t £ [t„ , i„+i ) , 

f(0,T) =0, 

for T £ [ti, T^+i), where a and b are the piecewise constant approximations 

a(i,r,x) |(t^^)g[t,,^t,,^j)x[rf,rj!+i) =a{tn,Te,x) = C^(a; + /o(i«)) X{tn,Te), 

K't^T,x) |(t,r)e[t,.,t„ + i)x[Tf,rf + i) =b{tn,Ti,x) ^ ^{x + fa{tn)) X{tn, Ti) . 

Thus, the extension above results in g{t, .) G for any time t e [0,t,nax]- 

Theorem 4.1. Letl :~ {0, . . . , A^}x{0, . . . , L— 1}, (g„ e)(n,e')€x be the numerical approximations produced 
by the (EFD) method. Also, we assume that the functions F , U , G along with their derivatives have 
polynomial growth. Also, we set 

(4.5) d{t,T,T,x) \ei.x + Ut))\{t,T)-\{t,T), 
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(4.2) 



(4.4) 



for x G M, t e [0,tinax] and r, r e [0,Tmax]- Then the computational error of the (EFD) method has the 
expansion 

:= E [g{g)] - E [^(1 )] = + Sz,,.„, + 0{{^tf + (Ar)^), 



(4.6) 

where 

(4.7) 

and 



Ed,, 



n=0 [ e=o 



2 ^n,e 



Ti=0 l t=0 £'=0 



2 full' 



n=0 



(4.8) 



£^x,,tim = E ^ 1 ^ [ ( J C^(^n+l) - F{gn,.) U{rn) ) ] 

+ E [ {fn+l - fn ) ^„+l,i,] 

+ X! ^ Wl)) ^(tn, n,Wn, ^n)))^n+l 

L-l 

X! ^ 1 X! ^ |^(^(i(t„+l,r^,T£/,f(<„+i,i„+i)) - d{tn,Tl,Tl-,"§{tn,tn))J$n+l.t.t 



ra=0 L f,f'=0 



with 



The two leading order terms i?r>,tau o,nd Eo^tim in the right hand side of (j4.6p are in a posteriori form 
and based on the discrete dualsTp^ £ R^+^ andTp^ G ]g(i-+2)x(i+2) ^/jj^/j (j,~g determined as follows. First, 
set 



i=l 



i=l 



for i = la,. ■ ■ , L — 1, and 



c„j(.t) := a(i„, r,-, x) Ai„ + b{t„, Tj,x) ■ AW„ 



for .T G M and j ~ 0, . . . , L — 1. Then, the first dual if is defined by the dual backward problem with final 
datum 



(4.9) 



and 



0, e-^ 

1, 



0,...,£a-l, 



(■a, ■ ■ ■ T L — \, 



L + l, 



(4.10) , = <^ 



?„+l,L + At„ F(f„ ,.) [/' (fn) ^n+l,L + l 

L-l 

^„+l,i + Atn ^"(f„,i,) 



^G{0,...,L-1}\{^4, 

£^ L, 
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0,...,4-l, ^' = 0,...,i + l, 



for n = N — I, . . . ,0. The second dual, , has final datum 

J G"(A,.c m A..cA^) A,,«,r (V), £,£' e {£a, . . . , L - 1}, £ ^ £' , 



(4.11) y'M,^' = { 



0, 
0, 

F'(^„,J G'(A,,Q(^))X,.«,r(^), 

F"(f„,J G(A,,«(f)), 

0, 



£ = £a,...,L-l, £' = L, 
£ = £a,...,L, £' 

£ = £aT ■ ■ ■ T L -\r I, £ = 0,...,^a — 1, 

£ = L, £' ^£a,...,L-l, 
£^L, (' = L, 

£ = L + 1, £' = £a,...,L + l, 



and solves the 



recursion 



f =' 



£,£'e{o,...,L-i}\{^„}, 



j",p=0 



2Ai„ X((5f„,j +4j(5„,c))'/'n+ijM- 

L-l 



+1 j'.i'+i 



(4.12) ^„,„, = <^ 



+ 2At„F(f„^jC/'(F„)PUi,.,.+i 

L-l 

+ Ai„ F(f„, J C/'(F„) f £ = 4, G {0, . . . , L - 1}\{£„} 

^G{0,...,L-l}\{f„}, £' = C, 



+ At„F'(g„^JC/(r„)^l+i^,^,^,,, ^ = L, G {0, . . . , L - 1}\{£„}, 

C^'f, €g{o,...,l-i}\{£„},/ = l, 



(4.13) ^„,,, 



(4.14) 



1+1,1.4-1 
L-l 

At„i^'(f„_jc/(f„)) X(<5,-,„ + c;,/f„,,j)p:,+i^ 



+ (At„)2F'(5„_J?7(r„)) +i^(5n,Jf^'(^«)<?L+i,L+i.L+i 

, (€,OG{(i,^n),(C,i)}, 



At. 
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(4.15) 

(4.16) 

and 



(4.17) 



+ (At„)2 (F'(f„,jC/(F„))'Pl+i,.+,,.+i, ^ = f = L, 
=PUi,i„i,+i + ^*n^'(fn,JC^(^")C+i,^+i,^+i, (AO e {(L,L+1),(L + 1,L)}, 



, ^ = L + 1, G {0, . . . , L - 1, L + 1}\{^„}, 



■1+1.L + 
L-1 



n-\-l,L + i,L 



+ Ai„ F(f „^ J C/'(F„) f ^ = i + 1, A = 

L 0n,i',e. € 6 {0, . . . , L - 1, L + 1}, A = L + 1. 

Proof. The proof is an application of Theorem 2.2 in [21j . To be able to split the time and maturity time 
discretization errors, introduce the semidiscretized fluxes a and b that, for ti < t < Te+i, are defined 
as a{t,T,x) = £_'^{x + fo{t)) X{t,Te), b{t,T,x) = £^{x + fo{t)) X{t,Te) and denote by g the corresponding 
semidiscrete in r solution. As a first step, replace the exact solution of (14.11) . g, by a finite dimen- 
sional approximation: a piecewise constant g*(t, ■), which is an Euler approximation with a much finer 
discretization, both in time t and maturity time r, than Thus, g* uses a time grid {tn)n=o much 
finer than {tn)n=oj ^^'^ ^ maturity time grid, (t£)^^q much finer than {Te)f^Q. Consequently, the num- 
ber of time steps satisfy P » N, M :» L, respectively, and At := maxo<m<p-itm+i — Un « Ai, 
At := maxo<m<j\/-i T'm+i — Tm « At. In the application of Theorem 2.2 in [21], include the t- 
discretization error terms a — a, & — 5 as well as the f-discretization terms a — a, 6 — 6 in the error 
expansion, following Lemmata 2.1-2.5 in [2], to obtain (|4.6ll4.8p for g replaced by the piecewise constant 
process 5* . For this purpose, observe that g can be also thought of as a piecewise constant function on 
the finer T-partition that defines 5*. The second step is to let A/, P 00 and Af, At — > 0, using 



max \g(t, ■) ~ g^(t, ■)\'^ + 



[g(t,T,„+i)-g(t,f„)]-[g.(t,f„+i)-g.(t,f„.)] 



O Af + (At> 



for t G [0, tmax] and m = 0, . . . , M — 2, along with similar estimates for the corresponding dual functions 
^, . . ., to control the higher order terms in the error expansion. The latter strong convergence estimates 
follow moving along the lines of the analysis of Section [H □ 

Remark 4.1. In the (EFD) method the t- discretization error of (|4.6p and (|4.7p can, by (14. 4p . (|4.5p . be 

expressed by 



^ J2 A<„ <^ ^ At, E [e{rn) ^„,J 



2t= \ = 1 \it„,Tt+l)-\{t„,Tl) 



(4.18) 



and the time discretization is 



L-l 

3 X! IE [^^(^n)^'n,«,r [A(i„,Tf+l) ■ A(t„,T£. + i) - A(i„,T£) ■ A(t„,Tf')] At^At,/ 



2 



(4.19) 



E [^^(5„+i,jC/(r„+i) - F(5„.Jt^(^n)] -f E [(r„+i - r„)^„+i_ 

L-l 

+ [(f (^n+l)A(t„+l,T,) -e'(^„)A(t„,T,))^„+l_, 

£=0 

1,-1 

+ 5 H E[(f (F„+i)A(f„+i,Tf) • A(i„+i,T£') 

- C^(^>i) A(t„,T£) • X{tn,Te')) 
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In Monte Carlo computations all the expected values in (|4.18p and (|4.19p are naturally approximated by 
sample averages. 

Remark 4.2. The analysis of the (EFE) method follows a similar line as the estimates of the (EFD) 
method. The difference lies in the t- discretization error, which by virtue of the orthogonality of both 
X — nX and X — HA to the subspace of piecewise constant functions S'a^j becomes second order accurate. 
Therefore, more careful expansions, including interpolation estimates, need to be carried out in order to 
capture the second order contributions from the t -discretization. 

5. Numerical experiments 

In this section we provide numerical evidence for the weak computational error (|2.12p of the numerical 
methods defined in Section [2] approximating the quantity of interest E [J^(/)] = E[t/((7)] described in 
(|1.4m.6p . In particular, we show results from numerical experiments with examples that have known 
exact solution which permit a straightforward derivation of an exact solution to compare with. The 
implementation uses double precision FORTRAN 77 and simulates the increments of the J independent 
Wiener processes by a double precision modification of the functions rani and gasdev proposed in [18] . 
The numerical quadrature approximation Q^'g) of A*(g) in (|2.9p is done via the use of Simpson's 
quadrature rule. For the particular case of the (EFD) method, the estimates for the computational 
error developed in Theorem 14.11 are compared with the exact computational error. The numerical results 
obtained are in agreement with the theory and the work to compute these estimates is small. 

5.1. Control of the statistical error. For M independent samples {Y{ujj)}j'^^ of a random variable 
y, with E [ |y 1^ ] < oo, define the sample average A{Y; M) and the sample standard deviation S{Y; M) 
of r by 



AI 



A{Y;M) = jjY.^iuj,) and SiY; M) = [AiY^; M) - {AiY; AI))^ y 

Let (T = ^E[|y — E[y]p] and consider the random variable 

Z„ = ^ (^(y;Af)-E[y]) 
with cumulative distribution function Fzj^^{x) = P{Zf,, < x), for x £ R. Let 

|3 " 



A 



1 



|y-E[y]r 



< CXD, 



then the Berry-Essecn theorem (cf. [10] p. 126), gives the following estimate in the central limit theorem 

sup|F,,,(x)-$(a;)| <-^X' 

for the rate of convergence of Fz^^^ to the distribution function, $, of a normal random variable with 
mean zero and variance one, i.e. 

*(-^) = ^ / 

Since in the examples below AI is sufficiently large, i.e. AI ^ 36 A^, the statistical error 

£:s(y;M) =E[y] -^(y;Af) 
satisfies, by the Berry-Esseen theorem, the following probability approximation 

p[[\£siY;AI)\<c, ^]) ^ 2$(co)-l. 
In practice choose some constant c„ > 1.65, so the normal distribution satisfies 

1 > 2$(cJ - 1 > 0.901 

and the event 

(5.1) |f,(y;M)|<E,(y;M)EEc, 

has probability close to one, which involves the additional step to approximate cr by S{Y;AI), cf. [TT] . 
Thus, in the computations Es{Y;M) is a good approximation of the statistical error £s{Y;AI). 
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For a given TOL > 0, the goal is to find M such that ¥ig{Y;M) < TOL. The algorithm described 
in |21| adaptively finds the number of realizations AI to compute the sample average AiY; M) as an 
approximation to ]E[K]. With large probability, depending on c^, the statistical error in the approximation 
is then bounded by TOL. For more details on the implementation of an adaptive algorithm to control 
the statistical error, see [21]. 

Remark 5.1 (Computational cost of the error estimates). The work to approximate E[^(g)] = E[X] 
within an accuracy TOL is O ^ ^^o^i ^ ; provided we use the Monte Carlo version of the EFD method as 
in (|2.11[) . It is therefore important to try to use both variance reduction techniques and adaptive methods 
to save computational effort. On the other hand, the work needed to compute sufficiently sharp error 
estimates as described in Theorem 14.11 is only ^^(TOL"^). The number of realizations needed to have a 

statistical error in the error bound much smaller than TOL is only ©(TOL^^) instead of the O ^ tol^^ ) 
realizations we need to compute an approximation of J-{g) using (|2.1ip . while the work to compute the 
error estimate for each realization is still 0{TOL~^), including the computation of the duals Tp and Tp . 
This surprising reduction of work for Tp and Tp is special for the HJM model studied here. For general 
SDEs the corresponding work would he ©(TOL^^) instead o/0(TOL^^). Thus, cheap and sharp error 
hounds are obtained hy the use of the a posteriori error estimates in Theorem 14. 11 Observe that if variance 
reduction techniques are applied to the approximation of ¥^[Q{g)], it is natural to try to use them also to 
reduce the variance in the error estimators. 

Remark 5.2 (Variance reduction techniques). The use of variance reduction techniques can decrease 
substantially the statistical errors. In particular the so called antithetic variates technique introduced in 
|12j reduces the variance in a sample estimator A{M; Y) by using another estimator A{M] Y) with the 
same expectation as the first one, but which is negatively correlated with the first. Then, the improved 
estimator is A{M; ^-i-^). Here, the choice ofY and Y relates to the Wiener process W and its reflection 
along the time axis, —W, which is also a Wiener process. If a realization of the Wiener process, W{-,U}j), 
yields, using one of the numerical discretizations p.5ll2.8p . a realization g(-,-,ujj) and —W{-,ijjj) yields 

'g{-,-,UJj) respectively, then we choose 

I ^ F{gN,L(^i)) G(A4,,Q(g(-,-,Wj)))+g^^^^(Ljj)+F(g„^(wj)) G (Ag;,,Q ,tjj ))) (wj ) 

as a better estimate. All the numerical results presented below use antithetic variates. In general, the use 
of control variates, see [5], can be also combined with other variance reduction methods. For example, 
the control variates technique is based on the knowledge of an estimator Y*-, positively correlated with Y , 
whose expected value IE[F*] is known and relatively close to the desired E[y], yielding Y — Yf, + E[l^] as 
an improved estimator. The estimates presented in this work do not preclude the use of control variates, 
and even though it is not applied here, it can be a valuable tool in practical computations. 

5.2. Numerical results. Now let us introduce some notation to be used later in the description of our 
numerical results. E^.^^ denotes the sample average approximating the r-discretization error (j4.7p and E^..^ 
denotes the sample average approximation to the i-discrctization error (j4.8l) . Beside this, denote by Eg the 
approximation (jS.ip to the statistical error Eg introduced in (|2.14p and by E^,^^ g the approximation (|5.ip 
to the statistical error in the estimation of the r-discretization error (|4.7|) by sample averages. Similarly, 
^tim s denotes the corresponding approximation to the statistical error in the estimation of the expected 
values in t-discretization error 



5.2.1. Ho-Lec model. The Ho-Lce model has ^(x) = a and \„{x) = 1 so A„ (x) = x and (|1.2p - (|1.3p takes 
the form 

df{t,T)^a'^{T-t)dt + adW{t), 0<t<T, 
^^■^^ /(0,r)=/„(r) 
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for T e [0, Tniax]- In this example the initial condition is /^(r) = r,-, — ^"J"^ + Jq '&{s)ds, where and a 
are real positive constants and d : M+ — M is a given function. Then, the exact solution of (|5.2p is 



^{T-tf+ d{s) ds + aW{t), 0<t<T, 



fit,r)=r 

which follows the normal distribution and therefore, yields bond prices which are log-normal distributed, 
allowing the use of Black and Scholes formulas for the pricing of call and put options on bonds. 

Setting Ta = ^max^ F{x) — 1 — X, G{x) = X, ^'(x) = X and U{x) = in (|1.5|) - (|1.6p . the functional to be 
computed has the form 



(5.3) 



E [T{f)] = E 



1 - 



/(s, s)ds 



/(imax, t) dr 



In the numerical experiments we choose r„ 0.05, a = 0.01, i?(.s) = -^e-". Then E [J'(/)] is a known 
function of imax and Tmax- The first experiment sets tmax = 1-0 and Tmax = 2.0, comparing the efficiency 
of the (EFD) and (EFE) methods. Table [523] shows the computational error for both methods and com- 
pares the a posteriori approximation of the error with the true computational error for the (EFD) method. 
Here, a confidence interval for the ratio between the error approximation and the exact computational 

-±^andi?-'^+'- 



B], with A 



When- 



error, S^., introduced in (|2.12p . is [A — B,A-, ^^j, ..^u^i _ uii^i _ |^-p 

ever we use the (EFD) method we call Sc.efd = and if we use the (EFE) method we call £c,efe = £c- 
Observe that the ratio A ± _B of the a posteriori approximation of the error over the computational error 
becomes closer and closer to one as we refine the time and maturity partitions, provided that the statis- 
tical error is small compared to the i-discretization error and the r-discretization error. In this example, 
the t-discretization gives the largest contribution to the computational error, and there is no practical 
advantage in the use of the (EFE) method. 



iseed = — 1 


(EFE) 


(EFD) 


N = L 




^ C.EFD 


[A- B,A + B] 


5 


-8.40 X 10'^ 


-8.25 X 10-* 


[0.97,0.97] 


10 


-4.16 X 10-^ 


-4.08 X 10"* 


[0.98,0.99] 


20 


-2.07 X 10--* 


-2.03 X 10"* 


[0.98,1.00] 



Table 15.2.11 Comparing the (EFD) and (EFE) methods in the Ho-Lee model approximating 
functional ([ESI) with M = 5000 and co = 1.65. 



e so 



5.2.2. Vasicek model. The Vasicek model has ^{x) = a and Ao(x) 

\(x) = ie—- (1-e-"-) 
and the forward rate equation (|1.2m.3p becomes 

df{t, r) = ^ (l - e-"(^"*)) e-O'ir-t) _^ ^ g-a(r-i) ^^r{^^)^ < t < T, 

/(0,r)==/„(r) 
for r G [0,Tijiax]- In this example the initial condition is 

h{r) = (ro - i) e-"^ + f - ^ (l - e""^)' , r G [t, w], 
where r^, a, a and § are given positive constants. The solution of (j5.4p is then 



(5.4) 



a e 
lo 



-a{t~~s) 



dW{s) 



1 



0<t<T, 



which is normally distributed and yields bond prices that are lognormal, as in the Ho-Lee model. 

Here we set = tmax = 0.3, Tmax = 6.0, and approximate again the functional defined in (j5.3p . In 
addition, we take = 0.03, a — 1.0, a — 0.01 and ^ — 0.05. Table [OH] displays the computational 
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errors for the (EFD) and (EFE) methods and compares the a posteriori approximatfon of the error with 
the true error for the (EFD) method. Observe that the ratfo A± B oi the a posterfori approximatfon of 
the error over the computational error becomes cfoser and cfoser to 1 as we refine the time and maturity 
partitions, provided that the statistical error is small compared to the t- and r-discretization error. 



iseed = — 1 


(EFE) 


(EFD) 


N 


^ C.EFE 


^C,EFD 


[A- B,A + B] 


5 


-2.30 X 10-'' 


-2.07 X 10-^ 


[1.92,1.95] 


10 


-2.05 X lO-'' 


-1.95 X 10-^ 


[1.03,1.05] 


20 


-1.06 X lO"'' 


-1.00 X 10-^ 


[0.99,1.02] 



Table [5221 Comparing the (EFD) and (EFE) methods in the Vasieek model approximating 
functional dO]) with M = 5000 and cq = 1.65. 



5.2.3. The Cox-Ingersoll-Ross (CIR) model. Consider the following (CIR) short rate model 
(5.5) r(t)=r, + / {d ~ ar{s)) ds+ [ a^/Hs) dW{s), t>0, 



where t?, a and a arc real constants. To connect the solution r(t) of (|5.5p to the diagonal value f{t,t) 
of the solution of an HJM problem, consider, first, the solution B = B{t; t) of the following Riccati 
differential equation (see [1]): 

dtB{t;T) :^^a^B^{t]T)+aB{t;T)-l, te[0,r], r > 0, 
B{t;t) =0, 

which has the form B{t] r) ~ tp{T — t) where 



sin 



m = -^ + ^ 7o ' 'I" . ' and 7o:=^y^^^T^- 



Provided ^(x) = (j-^max{x, 0} and \ — ip'{x)^ then A„(x) — ip'{x)ip{x) and the stochastic function 

f{t,T)^r{t)i''{T-t) + ^ijiT-t) 

solves (|1.2p - (|1.3p with the initial condition (t) = r,, ?/''(''') + '&iP{t). Taking into account that ip'{Q) ~ 1 
and '0(0) = 0, it follows that !{t,t) = r(t). 

Setting Ta = tmax, F(x) = e"^, G{x) = maxfe"-^ - K„,0}, *(.t) = x and U{x) = in (|1.5m.6p . the 
functional to compute in this example takes the form 

(5.6) E[n/)]=IE exp(^-^ /(s, s) ds ^ max jexp /(t„,ax, r) dr^ - 

In the numerical experiments we choose r„ = 0.15, a = 1.0, a = 0.1, = 0.05, <max = 5.0, Tmax = 8.0 
and Kg = 0.5. Table 15.2.31 shows the computational errors for the (EFD) and (EFE) methods and 
the ratio between the approximation of the computational error and the exact computational error for 
(EFD) method. There is no practical difference in this case between the (EFD) and the (EFE) method 
since the computational error is mainly t-discretization error and the r-discretization error is relatively 
unimportant. 



In order to have smooth coefficients in the HJM model (|1.2m.3p we approximate the function -y/maxja;, 0} 
in the diffusion term by a Lipschitz function globally defined in M (cf. [9] p. 252), 



v/max{2?70} w \J ^{x + \/x'^ + S) 



where 5 is a small positive constant. Observe that after this regularization the value of the functional 
E[J^(/)] depends on i5. In the computations 6 has been taken small enough to make this dependence 
negligible with respect to the size of the computational error. In this example we compute an accurate 
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iseed = — 1 


(EFE) 


(EFD) 


N = L 




^C,EFD 


[A-B,A + B] 


5 


1.23 X 10-^ 


1.21 X 10-^ 


[0.31,0.44] 


10 


5.83 X 10-^ 


5.39 X 10-^ 


[0.91,0.95] 


20 


2.76 X 10-^ 


2.79 X 10-^ 


[0.89,0.94] 



Table 15.2.31 Comparing the (EFD) and (EFE) methods in the (CIR) model approximating 
functional dHS]) with M = 2000 and cq = 1.65. 



numerical approximation of the exact E [J-{f) ] from (|5.6p . via the Feynman-Kac representation formula, 
using a numerical solution of the following backward PDF (cf. [201 P- 313), 

Vt + {'&-ar)vr + ^a'^rVrr-rv^Q, tG[0,imax], ?'e[0,rmax], 

with final datum u(imax,?') = (-^(r, ^max, Tmax) - Kg)~^ , where i?(r, i^ax, Tmax) denotes the (CIR) value 
for a bond with contracting time tmax, maturity time Tmax and short rate at tmax equal to r. We also use 
the boundary conditions 

vtit,0) + avrit,0)=0, i'(t, r,„ax) = 0, 

for t S [0,iinax]- The value of r^ax » ^ is taken sufficiently large so that the homogeneous Dirichlet 
boundary at r = Tmax has a negligible effect on the numerical approximation for v(0, 0.15) = E [J-'{f) ]■ 
The spatial discretization is a centered finite differences scheme and the time stepping is done by a 
diagonally implicit Runge Kutta method, namely the DIRK2 method, see [5]. Another way to estimate 
the exact solution with high accuracy is to use a formula based on the distribution (see |19| . pp. 
187-193 for details). 

5.2.4. A two-factor Gaussian model. A two-factor model has randomness introduced by two scalar in- 
dependent Wiener processes Wi, W2- In particular, for a two- factor Gaussian model we have £,{x) = 1, 
= '''I and Aq 2 (2;) = <y2e where cri, (72 and 02 are real positive constants. Thus (|1. 2^ - 1)1. 3p 



^0,1 V'*'/ 'Vq 2 

takes the form 



d/(t,r) 

(5.7) 



, ,2 -"2(--') 
2(^2)" 



(-1)^ {r-t)+ ^^^^' - - 1-e 



°2(T-t) 



dt 



CTi dWi[t) + a2e~'"'^^'^ *' dW2{t), 0<t<T, 



/(0,t)=/„(t) 



for T e [0,Tmax]- Here the initial condition is /(,(t) = fep + ^1 e^'^'^ where 6o: bi and k are real constants. 
Then, the exact solution of (|5.7p is normal distributed as in the Ho-Lee and Vasicek models, so explicit 
formulas are available for the pricing of put and call options with bonds as underlyings. 

In the numerical experiment we take ai ~ 0.02, 02 = 0.01, 02 = 0.5, and compute with the functional 
defined in (j5.6|) with strike K^^ = 0.5, tmax = 1 and Tmax = 3. For the initial condition we set 60 = 0.0759, 

= -0.0439 and k = 0.4454. Table [Qlil shows the computational errors for the (EFD) and (EFE) 
methods and the ratio between the approximation of the computational error and the exact computational 
error for method (EFD). 



iseed = — 1 


(EFE) 


(EFD) 


N 


^C.EFE 


^C^EFD 


[A- B,A + B] 


5 


-5.15 X 10-* 


-6.90 X 10-* 


[0.98,1.02] 


10 


-2.78 X 10-4 


-3.50 X 10-4 


[0.96,1.05] 



Table [1231 Comparing the (EFD) and (EFE) methods in the two-factor Gaussian model 
approximating functional (|5.6p with M — 40000 and cq = 1.65. 
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